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THEORY OF THE LATTICE BOLTZMANN METHOD: DISPERSION, DISSIPATION, 
ISOTROPY, GALILEAN INVARIANCE, AND STABILITY 

PIERRE LALLEMAND* AND LI-SHI LUO+ 


Abstract. The generalized hydrodynamics (the wave vector dependence of the transport coefficients) of 
a generalized lattice Boltzmann equation (LBE) is studied in detail. The generalized lattice Boltzmann equa- 
tion is constructed in moment space rather than in discrete velocity space. The generalized hydrodynamics 
of the model is obtained by solving the dispersion equation of the linearized LBE either analytically by using 
perturbation technique or numerically. The proposed LBE model has a maximum number of adjustable 
parameters for the given set of discrete velocities. Generalized hydrodynamics characterizes dispersion, dis- 
sipation (hyper- viscosities), anisotropy, and lack of Galilean invariance of the model, and can be applied to 
select the values of the adjustable parameters which optimize the properties of the model. The proposed 
generalized hydrodynamic analysis also provides some insights into stability and proper initial conditions for 
LBE simulations. The stability properties of some 2D LBE models are analyzed and compared with each 
other in the parameter space of the mean streaming velocity and the viscous relaxation time. The procedure 
described in this work can be applied to analyze other LBE models. As examples, LBE models with various 
interpolation schemes are analyzed. Numerical results on shear flow with an initially discontinuous veloc- 
ity profile (shock) with or without a constant streaming velocity are shown to demonstrate the dispersion 
effects in the LBE model; the results compare favorably with our theoretical analysis. We also show that 
whereas linear analysis of the LBE evolution operator is equivalent to Chapman-Enskog analysis in the long 
wave-length limit (wave vector k = 0), it can also provide results for large values of k. Such results are 
important for the stability and other hydrodynamic properties of the LBE method and cannot be obtained 
through Chapman-Enskog analysis. 

Key words, kinetic method, lattice Boltzmann equation, derivation of hydrodynamic equation, stability 
analysis, numerical artifacts of the LBE method 

Subject classification. Physical Sciences 

1. Introduction. The method of lattice Boltzmann equation (LBE) is an innovative numerical method 
based on kinetic theory to simulate various hydrodynamic systems [34, 5, 36]. Although the LBE method 
was developed only a decade ago, it has attracted significant attention recently [3, 6], especially in the area of 
complex fluids including multi-phase fluids [40, 41, 23, 32, 24, 25], suspensions in fluid [35], and visco-elastic 
fluids [12, 13]. The lattice Boltzmann equation was introduced to overcome some serious deficiencies of its 
historic predecessor: the lattice gas automata (LGA) [10, 46, 11]. The lattice Boltzmann equation circum- 
vents two major shortcomings of the lattice gas automata: intrinsic noise and limited values of transport 
coefficients, both due to the Boolean nature of the LGA method. However, despite the notable success of the 
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LBE method in simulating laminar [27, 31, 16, 17] and turbulent [45] flows, understanding of some impor- 
tant theoretical aspects of the LBE method, such as the stability of the LBE method, is still lacking. It was 
only very recently that the formal connections between the lattice Boltzmann equation and the continuous 
Boltzmann equation [19, 20, 1] and other kinetic schemes [28] were established. 

In this work we intend to study two important aspects of the LBE method which have not been sys- 
tematically studied yet: (a) the dispersion effects due to the presence of a lattice space; (b) conditions for 
stability. We first construct a LBE model in moment space based upon the generalized lattice Boltzmann 
equation due to d’Humieres [8]. The proposed model has a maximum number of adjustable parameters 
allowed by the freedom provided by a given discrete velocity set. These adjustable parameters are used to 
optimize the properties of the model through a systematic analysis of the generalized hydrodynamics of the 
model. Generalized hydrodynamics characterizes dispersion, dissipation (hyper- viscosities), anisotropy, lack 
of Galilean invariance, and instability of the LBE models in general. The proposed generalized hydrodynamic 
analysis enables us to improve the properties of the models in general. The analysis also provides us better 
insights into the conditions under which the LBE method is applicable and comparable to conventional CFD 
techniques. 

Furthermore, from a theoretical perspective, we would like to argue that our approach can circumvent 
the Chapman-Enskog analysis to obtain the macroscopic equations from the LBE models [8, 12, 13]. The 
essence of our argument is that the validity of the Chapman-Enskog analysis is entirely based upon the fact 
that there are two disparate spatial scales in real fluids: the kinetic (mean-free-path) and the hydrodynamic 
scale the ratio of which is the Knudsen number. When the LBE method is used to simulate hydrodynamic 
motion over a few lattice spacings, there is no such separation of the two scales. Therefore, the applicability 
of Chapman-Enskog analysis to the LBE models might become dubious. Under the circumstances, analyzing 
the generalized hydrodynamics of the model becomes not only appropriate but also necessary. 

It should also be pointed out that there exists previous work on the generalized hydrodynamics of the 
LG A models [33, 30, 15, 14, 7] and the LBE models [2], However, the previous work only provides analysis 
on non-hydrodynamic behavior of the models at finite wave-length, without addressing important issues 
such as the instability of the LBE method or providing insights as how to construct better models. In the 
present work, by using a model with as many adjustable parameters as possible, we analyze the generalized 
hydrodynamics of the model so that we can identify the causes of certain non-hydrodynamic behavior, such 
as anisotropy, and lack of Galilean invariance, and instability. Therefore, the analysis shows how to improve 
the model in a systematic and coherent fashion. 

This paper is organized as follows: Sec. 2 gives a brief introduction of the two-dimensional 9-velocity 
LBE model in discrete velocity space. Sec. 3 discusses the generalized LBE model in moment space. Sec. 4 
derives the linearized lattice Boltzmann equation from the generalized LBE model. Sec. 5 analyzes the 
hydrodynamic modes of the linearized evolution operator of the generalized LBE model, and the generalized 
hydrodynamics of the model. The dispersion, dissipation, isotropy, and Galilean invariance of the model are 
discussed. The eigenvalue problem of the linearized evolution operator is solved analytically and numerically. 
Sec. 6 analyzes the stability of the LBE model with BGK approximation, and compares with the stability of 
the LBE model presented in this paper. Sec. 7 discusses the correct initial conditions in the LBE simulations, 
and presents numerical tests of shear flows with discontinuities in the initial velocity profile. Sec. 8 provides 
a summary and concludes the paper. Two appendices provide additional analysis for variations of the LBE 
models. Appendix A analyzes a model with coupling between density p and velocity u , and Appendix B 
analyzes the LBE models with various interpolation schemes. 
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2. 2D 9- Velocity LBE Model. The guiding principle of the LBE models is to construct a dynamical 
system on a simple lattice of high symmetry (mostly square in 2D and cubic in 3D) involving a number of 
quantities which can be interpreted as the single particle distribution functions of fictitious particles on the 
links of the lattice. These quantities then evolve in a discrete time according to certain rules that are chosen 
to attain some desirable macroscopic behavior which emerges at scales large relative to the lattice spacing. 
One possible “desirable behavior” is that of a compressible thermal or athermal viscous fluid. (For simplicity 
of the analysis, we shall restrict our analysis to the athermal case in this work.) We shall demonstrate that 
the LBE models can satisfactorily mimic the fluid behavior to an extent that the models are indeed useful 
to simulate flows according to the similarity principle of fluid mechanics. For the sake of simplicity, we limit 
our discussions here in two-dimensional space. The extension to three-dimensional space is straightforward, 
albeit tedious. 

A particular two-dimensional LBE model considered in this work is the 9- velocity model. In this model, 
space is discretized into square lattice, and there are nine discrete velocities given by: 


(2.1) 


&01 


' ( 0 , 0 ) , 

< (cos[(a: — 1)7t/2], sin[(a — l)7r/2])c, 
(cos[(2a - 9)tt/4], sin[(2a - 9)7t /4])-\/2c, 


a = 0, 
a = 1-4, 
a = 5-8, 


where c = S x /S t is the unit of velocity, and S x and 6 t are the lattice constant of the lattice space and the 
unit of time (time step), respectively. From now on we shall use the units of 5 X = 1 and S t = 1 such that 
all the relevant quantities are dimensionless. The above discrete velocities correspond to the particle motion 
from a lattice node r,j to either itself, one of the 4 nearest neighbors (a = 1 -4), or one of the 4 next-nearest 
neighbors (a = 5-8). This model can be easily extended to include more discrete velocities and in space of 
higher dimensions, thus to include further distant neighbors where the particles move to in one time step. 
Nevertheless, “hopping” to a neighbor on the lattice induces inherent limitations in the discretization of 
velocity space. 

For the particular model discussed here, nine real numbers describe the medium at each node rj of a 
square lattice: 


{fa(rj)\a = 0, 1,..., 8}. 

The number f a can be considered as the distribution function of velocity e a at location r j (and at a particular 
time t). The set {/„} can be represented by a vector in R 9 which defines the state of the medium at each 
lattice node: 

(2-2) |/(r,)> = (/o, A, ...,/ 8 ) t . 

Once the vector \f(rj)) is given at a point rj in space, the state of the medium at this point is fully specified. 

The evolution of the medium occurs at discrete times t = ndt, (with <5 t = 1). The evolution consists of 
two steps: 

1. Motion to the relevant neighbors (modeling of advection); 

2. Redistribution of the {f a } at each nodes (modeling of collisions). 

These steps are described by the equation 

(2.3) fa(rj + e a ,t- 1-1) = f a (rj, t) + Q a (f) ■ 
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The above equation is the so-called Lattice Boltzmann equation (LBE). The lattice Boltzmann equation can 
be rewritten in a concise vector form: 

(2-4) \f( r j + e a ,t + 1)) = | /(rj, t )) -I- |A/) , 

where the following notations are adopted: 

(2.5a) \f{rj +e a ,t+ 1)) = (f 0 (rj +e 0 ,t+ 1), /i (rj +ei,t + 1), fs(r :j + e 8 ,f + 1)) T , 

(2.5b) |A/) = (fi 0 (/), fir (/),--., 0 8 (/)) t , 

so that | f(rj + e a ,t + 1)) is the vector of a state after advection, and |A/) is the vector of the changes in 
|/) due to collision Q. 

The advection is straightforward in the LBE models. The collisions represented by the operator may 
be rather complicated. However, f 1 must satisfy conservation laws and be compatible with the symmetry of 
the model (the underlying lattice space). This might simplify ft considerably. One simple collision model is 
the BGK model [4, 5, 36]: 

(2.6) ^ = ~[/a - /< eq) ] , 

where r is the relaxation time in unit of time step S t (which is set to be 1 here), and fa q ' is the equilibrium 
distribution function which satisfies the following conservation conditions for an athermal medium: 

(2.7a) d = E/« eq) =E/«> 

a oc 

(2.7b) pu = Y j e a f ( : q) = E e “/«’ 

a a 

where p and u are the (mass) density and the velocity of the medium at each lattice node, respectively. For 
the so-called 9-velocity BGK model, the equilibrium is usually taken as: 

(2.8) /£ eq) = w a p 1 + 3(e a • «) + |(e a • uf - ^ u 2 , 

where Wq = 4/9, uq, 2 , 3,4 = 1/9, and u’ 5 , 6 , 7,8 = 1/36. 

Some shortcomings of the BGK model are apparent. For instance, because the model relies on a single 
relaxation parameter r, the Prandtl number must be unity when the model is applied to thermal fluids, 
among other things. One way to overcome these shortcomings of the BGK LBE model [5, 36] is to use a 
generalized LBE model which nevertheless retains the simplicity and computational efficiency of the BGK 
LBE model. 

3. Moment Representation and Generalized 2D LBE. Given a set of b discrete velocities, 
{e a |a = 0, 1, ...,(&— 1)} with corresponding distribution functions, {f a \ot = 0, 1 ,...,(& — 1)}, one can 
construct a ^-dimensional vector space M 6 based upon the discrete velocity set, and this is usually the space 
mostly used in the previous discussion of the LBE models. One can also construct a space based upon the 
(velocity) moments of {/ a }. Obviously, there are b independent moments for the discrete velocity set. The 
reason in favor of using the moment-representation is somewhat obvious. It is well understood in the context 
of kinetic theory that various physical processes in fluids, such viscous transport, can be approximantly 
described by coupling or interaction among ‘modes’ (of the collision operator), and these modes are di- 
rectly related to the moments ( e.g ., the hydrodynamic modes are linear combinations of mass, and momenta 
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moments). Thus the moment-representation provides a convenient and effective means to incorporate the 
physics into the LBE models. Because the physical significance of the moments is obvious (hydrodynamic 
quantities and their fluxes, etc.), the relaxation parameters of the moments are directly related to the various 
transport coefficients. This mechanism allows us to control each mode independently. This also overcomes 
some obvious deficiencies of the usual BGK LBE model, such as a fixed Prandtl number, which is due to a 
single relaxation parameter of the model. 

For the 9-velocity LBE model, we choose following moments to represent the model: 


(3.1a) 

Ip) = (l, l, l, l, l, l, l, l, i) T , 

(3.1b) 

|e> = (-4, -1, -1, -1, -1,2, 2, 2, 2) 

(3.1c) 

\s) = (4, 2, 2, 2, 2, 1, 1, 1, 1) T , 

(3-ld) 

\jx) = (0, 1,0, -1,0, 1, -1, -1, 1) T , 

(3-le) 

\q x ) = (0, -2, 0, 2, 0, 1, -1, -1, 1) T , 

(3-lf) 

\jy) = (0, 0 , 1 , 0 , -1, 1, 1 , -1, -1) T , 

(3-lg) 

|9„> = (0, 0, -2, 0, 2, 1, 1, -1, -1) T , 

(3-lh) 

\Pxx) = (0, 1, -1, 1, -1, 0, 0, 0, 0) T , 

(3-li) 

|p*„> = (0, 0 , 0 , 0 , 0 , 1 , -1, 1 , -1) T . 


The above vectors are represented in the space ¥ = l 9 spanned by the discrete velocities {e a }, and they are 
mutually orthogonal to each other. These vectors are not normalized; this makes the algebraic expressions 
involving these vectors which follow simpler. Note that the above vectors have an explicit physical significance 
related to the moments of {f a } in discrete velocity space: | p) is the density mode; |e) is the energy mode; 
|e) is related to energy square; \j x ) and \j y ) correspond to the x- and y-component of momentum (mass 
flux); | q x ) and \q y ) correspond to the x- and y-component of energy flux; and \p xx ) and | p xy ) correspond to 
the diagonal and off-diagonal component of the stress tensor. The components of these vectors in discrete 
velocity space V = M 9 are constructed as follows: 


(3.2a) 

1 P)a 

= |e«|° = 1, 

(3.2b) 

| C ) a 

+ 

H 

CO 

+ 

o 

8 

V 

1 

II 

(3.2c) 

l £ )a 

= 4|e a |° - y(4,* + Ca,») + l ( e i 

(3. 2d) 

1 jx)a 


(3.2e) 

1 qx) ot 

= [— 5|e„|° + 3(e^ x + e 2 a<y )] e a , x , 

(3.2f) 

\jy)a 

— £a,y> 

(3-2g) 

\%)a 

= [-5|e a |° + 3{el tX + e 2 ^)] e a<y , 

(3.2h) 

| Pxx) a 

= e 2 - e 2 

^ a,x ot,y> 

(3.2i) 

1 Pxy) a 


Thus, 

(3.3a) 


p = (p\f) = ( f\p ) , 

(3.3b) 


e = (e|/> = </|e) , 

(3.3c) 


e = ( e l/) = (f\s), 
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(3.3d) 

(3.3e) 

(3.3f) 

(3-3g) 

(3.3h) 

(3-3i) 


jx — O'* I/) — (f\jx) 1 

Qx = ( q x \f ) = ( f\q x >, 

jy = (jy\f) = (f\jy) > 

tfy = 0,1/) = < f\Qy ), 

Pxx — (Pxx\f) — (f\P XX ) 1 
Pxy = {Pxy\f) = (f\Pxy) ■ 


Similar to {f a }, the set of the above moments can also be concisely represented by a vector: 
(3.4) |p) = (p, e, 6, j x , Qxi jy i Qyi Pxx i Pxy) ■ 

There obviously exists a transformation matrix M between |p) and |/) such that: 


(3.5a) k) = M|/), 

(3.5b) |/> = . 


In other words, the matrix M transforms a vector in the vector space V spanned by the discrete velocities 
into a vector in the vector space M = M 6 spanned by the moments of |/ 0 }. The transformation matrix M 
is explicitly given by: 


(3.6) 


/ (p\ ) 


/ 1 

1 

1 

1 

1 

1 

1 

1 

1 \ 

(e\ 


-4 

-1 

-1 

-1 

-1 

2 

2 

2 

2 

<e| 


4 

-2 

-2 

-2 

-2 

1 

1 

1 

1 

(jx | 


0 

1 

0 

-1 

0 

1 

-1 

-1 

1 

(Qx | 

= 

0 

-2 

0 

2 

0 

1 

-1 

-1 

1 

0,1 


0 

0 

1 

0 

-1 

1 

1 

-1 

-1 

0,1 


0 

0 

-2 

0 

2 

1 

1 

-1 

-1 

( Pxx | 


0 

1 

-1 

1 

-1 

0 

0 

0 

0 

\ (p.r/u 1 / 


l 0 

0 

0 

0 

0 

1 

-1 

1 

-1 / 


= (|p), |e), |e), | j x ), \q x ), \j y ), \ q y ), \p xx ), \p xy )) T ■ 


The rows of the transformation matrix M are organized in the order of the corresponding tensor, rather than 
in the order of the corresponding moment. The first three rows of M correspond to p, e, and s, which are 
scalars or zeroth-order tensors, and they are zeroth, second, and fourth order moment of {/„}, respectively. 
The next four rows correspond to j x , q x , j y , and q y , which are vectors or first-order tensors, and j x and j y are 
the first order moments, whereas q x and q y are the third order ones. The last two rows represent the stress 
tensor, which are second order moments and second order tensors. Again, this can easily be generalized to 
models using a larger discrete velocity set, and thus higher order moments, and in three-dimensional space. 
The main difficulty when using the LBE method to simulate a real isotropic fluid is how to systematically 
eliminate as much as possible the effects due to the symmetry of the underlying lattice. We shall proceed 
to analyze some simple (but non-trivial) hydrodynamic situations, and to make the flows as independent of 
the lattice symmetry as possible. 

Because the medium simulated by the model is athermal, the only conserved quantities in the system are 
density p and linear momentum j = (j x , j y ). Collisions do not change the conserved quantities. Therefore, 
in the moment space M, collisions have no effect on these three quantities. We should stress that the 
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conservation of energy is not considered here because the model is constructed to simulate an athermal 
medium. Moreover we find that the 9-velocity model is inadequate to simulate a thermal medium because 
it cannot have an isotropic Fourier law for the diffusion of the heat. Although the conserved moments are 
not affected by collisions, the non-conserved moments are affected by collisions, which in turn cause changes 
in the gradients or fluxes of the conserved moments which are higher order moments. In what follows the 
modeling of the changes of the non-conserved moments is described. 

Inspired by the kinetic theory for Maxwell molecules [26], we assume that the non-conserved moments 
relax linearly towards their equilibrium values that are functions of the conserved quantities. The relaxation 
equations for the non-conserved moments are prescribed as follows: 


(3.7a) 

e* = e — S 2 [e — e^ eq ^], 

(3.7b) 

e* = e - s 3 [s — £^ eq ^], 

(3.7c) 

Q* x =Qx- ■% [<lx - 4 eq) L 

(3.7d) 

<ly =%- «7 [q v - 4 eq) ], 

(3.7e) 

P*xx = Pxx - «8 [Pxx - 7>i e x q) ], 

(3.7f) 

Ply = Pxy ~ S 9 [Pxy ~ P^], 


where the quantities with and without superscript * are post-collision and pre-collision values, respectively. 
The equilibrium values of the non-conserved moments in the above equations can be chosen at will provided 
that the symmetry of the problem is respected. We choose 


(3.8a) 


(3.8b) 

(3.8c) 

(3.8d) 

(3.8e) 

(3.8f) 


e (eq ) = 


P. 


(e|e) 

1 


[«2 ( P\P ) P + 72 {{jx\jx)jl + {jy\jy)j 2 y)\ 

1 


r(eq) 


- 2 p + -72 {jx + jy), 

1 


e™ = 


[<*3 (p\p) P + 74 {(jx \jx)j x + {jy\jy)j 2 y)] 


7 (eq) _ 


(£|s) 

= -a 3 p+ -74 til + jy), 

{jx \jx) _ . _ 1 _ ■ 

(QxWx) 2 


0 (eq) 

Hy 


p (e q ) 

rXX 


_ tiy\jy) 
tiytiy) 

= 7i 


1 

Cl Jy — ~^ c l3y, 


( | A tix\jx)jl ~ tiy\jy)j 2 y ) = \lltil~il), 

\jr XX \Ir XX f 


(eq) = 

xy 




The values of the coefficients in the above equilibriums (ci, « 2 , 3 . and 71 , 2 , 3 , 4 ) will be determined in the next 
Section and summarized in Subsection 5.5. The choices of the above equilibriums are made based upon the 
inspection of the corresponding moments given by Eqs. (3.2), or the physical significance of these moments. 
Note that in principle q x and q y can include terms involving third order terms in terms of moment, such 
as jl and j x p xx [13], and e can include fourth order terms. Nevertheless, for the 9- velocity model, these 
terms of higher order are not considered because either they do not affect the hydrodynamics of the model 
significantly, or they lead to some highly anisotropic behavior which are undesirable for the LBE modeling 
of hydrodynamics. 
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Clearly LBE modeling of fluids is rather different from real molecular dynamics. Therefore it is not 
necessary to try and solve the mathematically difficult problem to create an inter-particle collision mechanism 
for the fictitious particles in the LBE models that would give the same eigenmodes of the collision operator 
in the continuous Boltzmann equation. However, what can be accomplished is that by carefully crafting a 
simple model with certain degrees of freedom, we can optimize large scale properties of the model in the 
sense that generalized hydrodynamic effects (deviations from hydrodynamics) are minimized. 

The values of the unknown parameters, c\, £* 2 , 3 , and 71 , 2 , 3,41 shall be determined by a study of the 
modes of the linearized collision operator with a periodic lattice of size N x x N y . 

It should be noted that in Eq. (3.8) the density p does not appear in the terms quadratic in j. This 
implies that the density fluctuation is decoupled from the momentum equation, similar to an incompressible 
LBE model with a modified equilibrium distribution function [18]: 


(3.9) 


/< eq) =^{ 


P + P 0 


3(e a -u) + -{e a - u) 2 



}' 


where the mean density po is usually set to be 1. The model corresponding to the equilibrium distribution 
function of Eq. (2.8) shall be analyzed in the Appendix A. 


4. Linearized LBE. We consider the particular situation where the state of the medium is a flow 
specified by uniform and steady density p (usually p = 1 so the uniform density may not appear in subsequent 
expressions) and velocity in Cartesian coordinates V = (V x , V y ) , with a small fluctuation superimposed: 

(4.1) I/) = |/ (0) > + \6f) 


where |/(°)) represents the uniform equilibrium state specified by uniform and steady density p and velocity 
y = (14 1 V y ), and | Sf) is the fluctuation. The linearized Boltzmann equation is: 

(4.2) | Sf(rj + e a , t + 1)) = \6f{r h t)) + Q (0) \6f(rj, t)) 


where is the linearized collision operator: 


(4.3) 


o(°) _ 

** “ w 0 


l/)=l/ <0) > 


*W({/< 0) }). 


In the moment space M, the linearized collision operator can be easily obtained by using Eqs. (3.7) and 
(3.8): 


(4.4) 


C pa — 


(Qp\qp) dA g a 


(Qa\Qa) dQ0 


le)=le (0) > 


where g a and |n Q ), a = 0, 1, .... b, are the moments defined by Eqs. (3.3) and the corresponding vectors 
in V = K 9 defined by Eqs. (3.1); A g a is the change of the moment due to collision given by Eqs. (3.7); 
|g) = |£>( 0 )} is the vector of all moments at the uniform equilibrium state [see Eq. (3.4) for the definition of 
|g)]. Obviously the linearized collision operator C depends on the uniform state specified by density p and 
velocity V = ( V x , V y ), upon which the perturbation |<5/) is superimposed. Specifically, for the 9-velocity 
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model, 

/ 0 0 0 0 0 0 000 
S2012/4: — s 2 0 .s' 27214/3 0 * 272^/3 0 0 0 

s 3 o; 3 /4 0 -s 3 S374I4/3 0 •s 3 7 -iK ,/3 0 0 0 

000 0 0 0 000 

( 4 . 5 ) C = 0 0 0 S5C1/2 -s 5 0 0 0 0 

000 0 0 0 000 

0 0 0 0 0 S7C1/2 — S7 0 0 

0 0 0 3s87iVi 0 — 3S871V), 0 — s 3 0 

V 0 0 0 3s 9 73^/2 0 3 s 9 7 3 14/2 0 0 -s 9 

The perturbation in the moments corresponding to \Sf) is |<$£), and |h^) = M|h/). The change of the 
perturbation due to collisions is linearly approximated by |A^) = C|(5p) in the moment space M spanned by 
{|<? a )|a = 0, 1, . . . , (6 — 1)}. This change of state in discrete velocity space ¥ is |A/) = M“ 1 C|(5^). Therefore 
the Eq. (4.2) becomes 

(4.6) | Sf(rj +e a ,t + 1)) = \8f(rj, t)) + M^CM^/^-, t )) . 

In Fourier space, the above equation becomes: 

(4.7) A| 8f(k, t + 1)) = [I + M _1 CM] | Sf(k, t )) , 

where A is advection operator represented by the following diagonal matrix in discrete velocity space ¥ = E 9 : 

(4.8) A a0 = exp(« e a ■ k)S a0 , 

where 6 a0 is the Kronecker 8. It should be noted that for a mode of wave number k = ( k x , k y ) in Cartesian 
coordinates, the advection operator A in the above equation can be written as follows: 

(4.9) A = diag(l, p, q, 1/p, 1/q, pq, q/p, 1/pq, p/q ) , 
where 

(4.10) p = e ik * , q = e iky . 

The advection can be decomposed into two parts, along two orthogonal directions, such as a:-axis and y-axis 
in Cartesian coordinates: 

A (k x ) = A {k x , k y = 0) = diag(l, p, 1, 1/p, 1, p, 1/p, 1/p, p ) , 

A (k y ) = A {k x = 0, k y ) = diag(l, 1, q, 1, 1/q, q, q, 1/q, 1/q) . 
and A (k x ) and A (k y ) commute with each other: 

A = A(k x )A(k y ) = /K(k y )A(k x ) , 

i.e., the advection operation can be applied along .-r-direction first, and then along .(/-direction, or vice versa. 
The linearized evolution equation (4.7) can be further written in a concise form: 

(4.11) \5f(k,t + l)) = L\5f(k,t)), 
where 

(4.12) L = A -1 [l + M -1 CM] , 
is the linearized evolution operator. 
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5. Modes of Linearized LBE. 


5.1. Hydrodynamic modes and transport coefficients. The evolution equation (4.6) is a difference 
equation which has a general solution: 

(5.1) \G( r j, t = l)) = K™ Ky Z l \Go) , 

where m and n are indices for space (r j = mx + ny), and x and y are units vectors along the .r-axis and 
y-axis, respectively; |Go) is the initial state. We can consider the particular case of a periodic system such 
that the spatial dependence of the above general solution can be chosen as 

(5.2) | Sf) = exp (—ik ■ rj + zt)\G(rj, t )) . 

By substituting Eqs. (5.1) and (5.2) into the linearized LBE (4.11), we obtain the following equation: 

(5-3) z\G 0 ) = L|G 0 ) , 

The above equation leads to the dispersion relation between 2 and fc: 

(5.4) det[L - z\] = 0 , 

which determines the transport behaviors of various modes depending on the wave vector k. The solution of 
the above eigenvalue problem of the linearized evolution operator L provides not only the dispersion relation, 
but also the solution of the initial value problem of Eq. (4.11): 

b 

I Sf(k, t)) = L*| <m 0)) = £ zi\z a )(z a \Sf(k, 0)) , 

a=l 

where | z a ) is the eigenvector of L with eigenvalues z a in discrete velocity space V. 

The eigenvalue problem of Eq. (5.4) cannot be solved analytically in general, except for some very special 
cases. Nevertheless, it can be easily solved numerically using various packages for linear algebra, such as 
LAPACK. For small fc, it can be solved by a series expansion in k. The only part in L which has fc-dependence 
is the advection operator A. Therefore, we can expand A -1 in L: 

(5.5) K = A" 1 = K<°> + K«(fc) + K < 2 >(fc 2 ) + • • • + K<")(fc n ) + • • • , 
where KW depends on k n : 

(5.6) K ^ = ±(-ik-e a ) n 6 a0 . 

When fc = 0, the eigenvalue problem of the ( b x 6)-matrix iJ 0 - 1 = (I + M -1 CM) can be solved analytically. 
There exists an eigenvalue of 1 with three-fold degeneracy, which corresponds to three hydrodynamic (con- 
served) modes in the system: one transverse (shear) and two longitudinal (sound) modes. It is interesting 
to note that when fc = (tt, 0) or fc = (0, n), L also has an eigenvalue of —1, which corresponds to the 
checkerboard mode, i.e., it is a conserved mode of L 2 . Being a. neutral mode as far as stability is concerned, 
it will be necessary to study how it is affected by a mean velocity V. Thus we shall have to analyze the 
model for fc ranging from 0 to tt, which the standard Chapman-Enskog analysis cannot do. 

The hydrodynamic modes at fc = 0 are: 

(5.7a) |0 t) = cos 0\j x ) - sm6\j y ) = \j T ), 

(5.7b) |p±) = |p) ± (cos% x ) +sin0|jg,)) = | p) ± | j L ), 
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where 6 is the polar angle of wave vector k. For finite fc, the behavior of these hydrodynamic modes depends 
upon k. In two-dimensional space, these linearized hydrodynamic modes behave as follows [29]: 

(5.8a) | grit)) = Zt\6t(0)) = exp [—ik(gV cos</>)f] exp(— vk 2 t)\gT(0)), 

(5.8b) \6±{t)} =-4|e±(0)) =exp[±ik(c s ±gV cos (j>)t]exp[-(v/2 + T])k 2 t\\g±(0)}, 

where v and g are the shear and bulk viscosity, respectively; the coefficient g indicates whether system is 
Galilean invariant (that g = 1 implies Galilean invariance); c s is the sound speed; V is the magnitude of the 
uniform streaming velocity of the system V = (V x , V y ); and <j> is angle between the streaming velocity V and 
the wave vector k. The Galilean-coefficient g(k) is similar to the ^-factor in the FHP lattice gas automata 
[10, 46, 11], which also determines the Galilean invariance of the system. The transport coefficients and the 
Galilean-coefficient are related to the eigenvalues of L as the following: 

(5.9a) v{k) = — y^-Re[ln zr(k)\, 

K 

(5.9b) <jr(fe)V cos <j> = — yIm[lnzT(fc)], 

k 

(5.9c) ^(fc) + v(k) = --^Re[ln z±(fe)], 

(5.9d) c a (k) ± g(k)V cos (j> = =F^Im[ln z±(fc)], 

A t 

where Zr(k) and z±(k) are the eigenvalues corresponding to the hydrodynamic modes of the linearized 
evolution operator L. Since the transport coefficients can be obtained through a perturbation analysis, we 
shall use the following series expansion in k: 


(5.10a) 

v(k) = vq - v\ k 2 + . . 

. + (-1 ) n v n k 2n + . 

(5.10b) 

rj(k) = go — r)\k 2 + .. 

. + (-1 ) n g n k 2n + . 

(5.10c) 

c t (k) = C 0 - Ci k 2 + . 

.. + (-1 ) n C n k 2n + 

(5.10d) 

#(&) = 9o — 9ik 2 + ■ ■ 

. + (-1 ) n g n k 2n + .. 


It should be noted that, in the usual Chapman-Enskog analysis of LBE models, one only obtains the values of 
the transport coefficients at k = 0. As we shall demonstrate later, higher order corrections to the transport 
coefficients (i.e., hyper- viscosities) are important to the LBE hydrodynamics, especially for spatial scales of 
a few lattice spacings. 

One possible method to solve the dispersion relation det[L — z\] = 0 is to apply the Gaussian elimination 
technique using l/s* as small parameters for the non-conserved modes (the kinetic modes). Starting from a 
9 x 9 (b x b in general) determinant, we obtain a 3 x 3 determinant for the 3 conserved modes. The elements 
of this new determinant are computed as series of 1 / Si and k with the necessary numbers of terms to achieve 
a given accuracy when computing the roots of the dispersion equation. 

It should be mentioned that the interest of the present technique is that it provides a very simple means 
to analyze models with various streaming and collision rules with as many adjustable parameters as possible 
to be determined later when trying to satisfy either the stability criteria or physical requirements to model 
various hydrodynamic systems. Free parameters are the equilibrium coefficients in Eqs. (3.8): Ci, cq, and 
7 i\ and relaxation rates Sj. 


5.2. Case with no streaming velocity ( V = 0). We first consider the case in which the streaming 
velocity V = 0. To the first order in k, we obtain two solutions of Im(ln z±) = =F ikc s with 


(5.11) 



li 



These are the sound modes supported by the medium. At the next order, we obtain modes with Re(ln zt) = 
—vq k 2 . To enforce isotropy we need to have 


(5.12) 


= (ci +4) 

s 9 2 2y (2 — ci) 


such that the ^-dependence in i/ 0 vanishes, 


(5.13) 


vo 


(2 ~ci) 1\ 

12 Vs 8 2 ) ' 


which can be interpreted as the shear viscosity of the medium in the limit k = 0 (measured in basic units of 
space and time). For the sound modes, we also find an attenuation rate Re(ln z-t) = — (i/ 0 /2 + rjo)k 2 where 
(z/ 0 /2 + rjo) is the longitudinal kinematic viscosity in a two-dimensional system. The bulk viscosity of the 
model at long wave length limit k = 0 is: 


(5.14) 


(ci + 10 - 12 C 2) / 1 1\ 

24 Vs 2 2 ) ‘ 


The positivity of the transport coefficients leads to the bounds on the adjustable parameters: 


(5.15a) -16 <£* 2 , 

(5.15b) —4 <ci < 2, 


and the bounds on the following relaxation parameters: 


(5.16a) 0 <s 2 < 2, 

(5.16b) 0 <s 8 < 2. 


The bounds for a 2 and C\ will be further narrowed in the following analysis. Based upon the above results 
of i/o, r? o, and c s , it is clear that the model is isotropic at rest (i.e., the streaming velocity V = 0) and in 
the limit of k = 0. The Galilean-coefficient g cannot be determined when the streaming velocity V = 0. 
Therefore, the case of a finite streaming velocity V is considered next. 

5.3. Case with a constant streaming velocity V. As indicated by Eqs. (5.9), to the first order 
in k, the three hydrodynamic roots of the dispersion equation (zt and z±) give the phase gV cos <j> and the 
sound speed c s . In order to make the root of the transverse mode (z-r) to have a correct phase corresponding 
to the streaming velocity V, as expected for a model satisfying Galilean invariance, i.e., go = 1, we must set 

2 

(5-17) 7i=73=g. 

If we further set 


(5.18) 72 = 18, 

then we obtain the roots of the sound modes (^±) which lead to the sound speed 

(5.19) C s = V cos (j> ± \/c 2 - \-V 2 cos 2 <j), 

where V cos </> = V ■ k, and k is the unit vector parallel to k. This clearly shows that the system obeys 
Galilean invariance only up to first order in V. One way to correct this defect is to allow for compressibility 
effects in the equilibrium properties, as shown in Appendix A. The dispersion of sound can be computed 
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either analytically, by carrying out the perturbation expansion in fc or numerically, by solving the eigenvalue 
problem for any value of k. The dispersion of sound is important when studying the nonlinear acoustic 
properties of the medium. 

Second, the attenuation of transverse wave depends not only on V but also on the direction of the wave 
vector k. In order to eliminate the anisotropy in the V-dependence of the shear wave attenuation, we must 
choose: 


(5.20) 


a = - 2 . 


With the above choice of C \ , the shear viscosity in the limit of k = 0 is given by: 

i/ 0 = [s 2 (2 - s 8 )(c 2 + (1 - 3c 2 )G 2 cos 2 4> ) + 3(2(s 8 - s 2 ) 

(5.21) +s 8 (s 2 - 2) cos 2 tp)V 4 cos 2 (/>]/[6s 2 ss(I /2 cos 2 </> + c 2 )] . 

Similarly, from the attenuation of acoustic waves, one obtains the bulk viscosity (in the limit of k = 0) which 
has a complicated dependence on the streaming velocity V : 

rio = {V cos 4>\/V 2 cos 2 tp + <? a [\2V 2 ((s 2 - s 8 ) + s 2 (s 8 - 2) cos 2 </>) 

+(2s 2 — 3s 2 s 8 + 4s 8 )(1 — 3c 2 )] 

(5.22) +3V 4 cos 2 (/>[cos 2 </>(2s 8 + 3s 2 s 8 - 8s 2 ) + 6 (s 2 - s 8 )] 

+2V 2 cos 2 <f>[6(s 2 s 8 -s 2 - S 8 )c 2 + Sg(2 - s 2 )] 

+c 2 [6V 2 (s 2 - s 8 ) + s 8 ( 2 - s 2 )( 2 - 3c 2 )]}/{12s 2 S8(^ 2 cos 2 </> + c 2 )} . 


It is obvious that the streaming velocity V has a second order effect on i/ 0 , and a first order effect on i) 0 . A 
careful inspection of the above result of r ) 0 indicates that the first order effect of V on t ] 0 can be eliminated 
by setting c 2 = 1/3 (or equivalently a 2 = -8). Furthermore, the second order effect of V on the sound 
speed and the longitudinal attenuation can also be eliminated by using a slightly more complicated model 
with thirteen velocities, as noted by a previous work [37]. 

In summary, although all the transport coefficients are isotropic in the limit fc = 0, some undesirable 
features of the LBE models can be clearly observed at the second order in k when the streaming velocity 
V has a finite magnitude. First, the acoustic wave propagation is not Galilean invariant. Second, both 
the shear and the bulk viscosity depend on V. Nevertheless, these effects are of the second order in V, 
and can be improved to higher order in both k and V by incorporating compressibility into the equilibrium 
properties of the moments (see Appendix A) or using models with a larger velocity set. 


5.4. Third order result. The analysis in the previous subsections shows that isotropy for the hydro- 
dynamic modes of the dispersion equation can be attained to the first and second order in k by carefully 
adjusting the parameters in the model. In the situation with a uniform streaming velocity V parallel to fc, 
we find that the third order term in fc for the shear mode is anisotropic, i.e., 


1 


3-s 2 


3s 8 


«i = -hr3-^ + « + o- — + j v cos ^ 


l 


«8 


(5.23) 


+ U-- + — (cos 4 0-cos 2 0 + f). 
3 s 8 s 5 \s 8 )\ \ 3 J 


The anisotropic term in gi (depending on cos 9) can be eliminated if we choose 

.(2 -ss) 


(5.24) 


■s 5 = 3 


(3 - 8s) 
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As indicated by Eq. (5.13), parameter s$ is usually chosen close to 2 from below in order to obtain a small 
shear viscosity (consequently large Reynolds number). Therefore, the preceding expression yields a small 
value for s 5 . This would lead to an undesirable consequence: Mode | q x ) relaxed with the relaxation parameter 
$5 would become a quasi-conserved mode leading to some sort of visco-elastic effect [13]. Therefore we usually 
choose to have large S5 such that the advection coefficient of transverse waves has an angular dependence 
for non-zero k in third order in k. That is, the physical conservation laws are preserved at the expense of 
the isotropy of the dispersion in third order (and all higher order) in k. 

It should be noted that the value of g has effects on the Reynolds number because the time t needs to 
be rescaled as gt. 

5.5. Optimization of the model and connection to the BGK LBE model. Among seven ad- 
justable parameters (ci, a*, and 7*) in the equilibrium values of the moments in the model [see Eqs. (3.8)], 
so far only five of these parameters have been fixed by enforcing the model to satisfy certain basic physics 
as shown in the previous analysis: c\ = —2, a 2 = —8, 71 = 73 = 2/3, and 72 = 18. These parameter values 
are the optimal choice in the sense that they yield the desirable properties (isotropy, Galilean invariance, 
etc.) to the highest order possible in wave vector k. It should be stressed that the constraints imposed by 
isotropy and Galilean invariance are beyond the conservation constraints — models with only conservation 
constraints would not necessarily be isotropic and Galilean invariant in general, as observed in some newly 
proposed LBE model for non-ideal gases [44, 43, 32]. Two other parameters, <*3 and 74, remain adjustable. 
In addition, there are six relaxation parameters s t in the model as opposed to one in the LBE BGK model. 
Two of them, s 2 and s$ determine the bulk and the shear viscosity, respectively. Also, because Ci = -2, 
therefore S9 = -s$ [see Eq. (5.12)]. The remaining three relaxation parameters, S3, S5, and Sj can be adjusted 
without having any effects on the transport coefficients in the order of k 2 . However, they do have effects in 
higher order terms. Therefore, one can keep values of these three relaxation parameters only slightly larger 
than 1 (no severe over-relaxation effects are produced by these modes) such that the corresponding kinetic 
modes are well separated from those modes more directly affecting hydrodynamic transport. 

It is interesting to note that the present model degenerates to the BGK LBE model [5, 36] if we use a 
single relaxation parameter for all the modes, i.e., s* = 1/r, and choose 

(5.25a) 0:3 = 4, 

(5.25b) 74 = -18. 

Therefore, in the BGK LBE model, all the modes relax with exactly the same relaxation parameter so there 
is no separation in time scales among the kinetic modes. This may severely affect the dynamics and the 
stability of the system, due to the coupling among these modes. 

6. Local Stability Analysis. The stability of the LBE method has not been well understood, although 
there exists some preliminary work [42, 47]. However, previous work does not provide much theoretical insight 
into either the causes or the remedies for the instability of the LBE method. In the following analysis, a 
systematic procedure which identifies some causes of instability is discussed and illustrated by some examples. 

Our stability analysis relies on the eigenvalue problem for the linearized evolution operator L, the disper- 
sion equation. For large values of k, one could in principle analyze the dispersion equation to higher order 
by perturbation expansion. In practice, it is more efficient to compute the roots of the dispersion equa- 
tion numerically. We shall try to identify the conditions under which one of the modes becomes unstable: 
instability occurs when Re(ln2 a ) < 0. 
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Table 6.1 

Special properties of the dispersion relation when the wave vector k is of some special values. 


k 

dispersion equation 

conditions 

(0, 0) 

[z - l] 3 = 0 
[Z ~ (1 - *2)] = 0 
[z - (1 - S 3 )] = 0 
[Z - (1 - s 5 )] 2 = 0 

[z ~ (1 - * 8 )] 2 = 0 

S7 = s 5 

(±1, 0)7 r 
or 

(0, ±1)t r 

[z + 1] = 0 

[z + (1 — s 5 )] = 0 or [z + (1 - s 7 )] = 0 
[z + (1 - s 8 )] = 0 or [z + (1 - s 9 )] = 0 
[ z 2 - | s 5 z + s 5 - 1] = 0 
[ 2 4 + |(s 3 - 2 s 2 )z 3 

+ |{s 2 (s 8 — 4s 3 ) — 6s 3 s 8 + 9(s 2 + S3 + s$ — 2)}z 2 
+ §(s 8 - 1X02(03 _ 2) + S 3 )z 
+ (1 - s 2 )(l - «3)(1 - $8)] = 0 


(±1,±1)tt 

[ z -( 1 -s 8 )] 2 = 0 
[ Z 2 - \s$z + s 5 - l] 2 = 0 
[z 3 + |(lls 2 - 3s 3 - 9 )z 2 
+ |{3(4s 3 - 3) - s 2 (s 3 + 2 )}z 
+ (1 - 0 2 )(1 - s 3 )] = 0 



We have noticed some interesting qualitative properties of the dispersion for the 9-velocity model when 
wave vector k is parallel to certain special directions with respect to the lattice line. These properties are 
listed in Table 6. These qualitative behaviors of the dispersion equation already demonstrate the strong 
anisotropy of the dispersion relations dictated by the lattice symmetry. 

To exhibit the complex behavior of the dispersion equation, we compute the roots of the dispersion 
equation with a given set of parameters. Figures l(a.) and 1(b) show the real and imaginary parts of the 
logarithm of the eigenvalues as functions of k, respectively. Figs. 1 clearly exhibit the coalescence and 
branching of the roots. This suggests a complicated interplay between the modes of collision operator 
affecting the stability of the model. The asymmetric feature of these curves is due to the presence of a 
constant streaming. 

The growth rate of a mode | z a ), Re(lnz a ), depends on all the adjustable parameters: the relaxation 
parameters, the streaming velocity V , and the wave vector k. To illustrate this dependence, we consider the 
BGK LBE model with 1/r = 1.99. Figure 2 shows the growth rate for the most unstable mode as a function 
of streaming velocity V and wave vector k. For each V, we let k be parallel to V, with a polar angle 0 
with respect to the .T-axis. Then we search for the most unstable mode in the interval 0 < k < tt. For the 
9-velocity BGK LBE model, the unstable mode starts to appear above V ~ 0.07. Figure 2 shows the strong 
anisotropy of the unstable mode: the growth rate significantly depends on the direction of fc, and the critical 
value of k at which the unstable mode starts to appear is also strongly anisotropic. We also compute the 
growth rate for the most unstable mode with V perpendicular to fe, and find that the stability of the model 
is generally qualitatively the same as when V is parallel to fe, but is slightly more stable. Generally we 
find that the transverse mode is more stable than longitudinal modes. In many instances we have observed 
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Fig. 1. The logarithmic eigenvalues of the 9-velocity model. The values of the parameters are: 0.2 = —8, 03 = 4, 
ci = —2, 71 = 73 = 2/3, 72 = 18, and 74 = —18. The relaxation parameters are: S 2 = 1.64, S3 = 1.54, S5 = S7 = 1.9, and 
ss = sg = 1.99. The streaming velocity V is parallel to k with V = 0.2, and k is along x-axis. (a) Re(ln^ a ), and (b) Im(lnz a ). 



V 

Fig. 2. The growth rate of the most unstable mode for the BGK LBE model — ln2 a vs. the streaming velocity magnitude 
V . The relaxation parameter ss = 1/r = 1.99. The wave vector k is set parallel to the streaming velocity V. For each value 
of V with a polar angle 6 with respect to x-axis, the growth rate is computed in the interval 0 < k < w in k-space. Each curve 
corresponds to the growth rate of the most unstable mode with a given V , and k parallel to V with the polar angle 0 with 
respect to x-axis. 


that sound waves propagating in the direction of the mean flow velocity V can be quite unstable. This 
instability may be reduced by making the first order V-dependent term in the attenuation of the sound 
waves [f?o in Eq. (5.22)] equal to 0 by choosing c 2 s = 1/3, as indicated in the previous section. It should be 
noted that when the growth rate is infinitesimal, it takes extremely long time for the instability to develop in 
simulations. Because the unstable modes we have observed have a large wave vector k (small spatial scale), 
therefore, as a practical means to reduce the effect of instabilities in LBE simulations, some kind of spatial 
or temporal filtering technique may be used in the LBE schemes to reduce small scale fluctuations and thus 
to limit the development of instabilities. 

It should be pointed out that we do not discuss here the influence of boundary conditions which may 
completely change the stability behavior of the model through either large scale genuine hydrodynamic 
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s 8 =l/r 

Fig. 3. The stability of the generalized LBE model vs. the BGK LBE model in the parameter space of V and s g = 1/r. 
The line with symbol □ and X are results for the BGK LBE model and the model proposed in this work, respectively . The 
region under a curve is the stable region in the parameter space of V and s% = 1/r. Note that the stability of the BGK LBE 
model starts to deteriorate after s$ > 1.92, whereas the stability of the proposed generalized LBE model remains virtually intact. 



behavior or local excitation of Knudsen modes. 

As previously indicated, the adjustable parameters in our model can be used to alter the properties of 
the model. The stability of the BGK LBE model and our model is compared in Fig. 3. In this case we 
choose the adjustable parameters in our model to be the same as the BGK LBE model, but maintain the 
freedom of different modes to relax with different relaxation parameters s*. Fig. 2 shows that for each given 
value of V, there exists a maximum value of s$ = 1/r (which determines the shear viscosity) below which 
there is no unstable mode. The values of other relaxation parameters used in our model are: s 2 = 1.63, 
S 3 = 1.14, s 5 = s 7 = 1.92, and Sg = Ss = 1/r. Fig. 3 clearly shows that our model is more stable than 
the BGK LBE model in the interval 1.9 < Sg = 1/r < 1.99. Therefore, we can conclude that by carefully 
separating the kinetic modes with different relaxation rates, we can indeed improve the stability of the LBE 
model significantly. 

7. Numerical Simulations of Shear Flow Decay. To illustrate the dispersion effects on the shear 
viscosity in hydrodynamic simulations using the LBE method, we conduct a series of numerical simulations 
of the shear flow decay with different initial velocity profiles. The numerical implementation of the model is 
discussed next. 

7.1. Numerical implementation and initial conditions. The evolution of the model still consists 
in two steps: advection and collision. The advection is executed in discrete velocity space, namely to 
\f(x, t)), but not to the moments | q(x, t)). However, the collision is executed in moment space. Therefore, 
the evolution involves transformation between discrete velocity space V and moment space M, similar to 
Fourier transform in the spectral or Galerkin methods. The evolution equation of the model is: 

(7.1) | f(x + e a 8 t , t + S t )) = \f(x, t)) + M _1 S [|p(cc, t)) - |p (eq) }] , 
where S is the diagonal relaxation matrix: 

(7.2) S = diag(0, -s 2 , -S 3 , 0, -s 5 , 0, -s 7 , -s 8 , -s 9 ) . 
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In simulations using the LBE method, the initial conditions provided are usually specified by velocity 
and pressure (density) fields. Often the initial condition of f a is set to its equilibrium value corresponding 
to the given flow fields, with a constant density if the initial pressure field is not specified. The initial 
conditions of f a can include the first order effect fa ] ■ The first order effect in moment space is obtained 
through Eq. (7.1): 

(7.3) \ e W) = S- 1 MD\fW), 
where D is a diagonal differential operator: 

(7.4) D a fj = 8 a 0 C a • V . 

Eq. (7.3) is similar to Chapman-Enskog analysis of fa ' ■ 

For the shear flow, only the initial velocity profile is given. The density mode is set to be uniform 
initially. The remaining modes are initialized as the following: 


(7.5a) 

P = 1> 

(7.5b) 

e = -2 + 3 (u 2 x + u y ). 

(7.5c) 

€ = l- 3(u x + u 2 y ), 

(7.5d) 

Qx — 

(7.5e) 

( 

53 

1 

II 

(7.5f) 

Pxx — ( u x u y ) g s iPx^x 9 y Uy), 

(7-5g) 

Pxy — 'U'x^y ^ i,^y^x ^x^y)‘ 

oS$ 


The terms in p xx and p xy involving derivatives of the velocity field take into account of viscous effect in 
the initial conditions. These terms are obtained through Eq. (7.3). The first order terms in turn induce 
second order contributions (with respect to space derivatives) which are not included here. This leads to 
weak transients of short duration if there is separation of time scales (2 - s 8 ) <C (2 — s 5 ). 

Our first test is the decay of a sinusoidal wave in a periodic system for various values of k. The numerical 
and theoretical results agree with each other extremely well and confirm the fc-dependence of g and v. The 
agreement indicates that our local analysis is indeed sufficiently accurate in this case. 

The next case considered is more interesting and revealing because the initial velocity contains shocks. 
Consider a periodic domain of size N x x N y = 84 x 4. At time t = 0, we take a shear wave u y (x, 0) of 
rectangular shape (discontinuities in u y at x = N x / 4 and x = 3N x /4): 

u y (x, 0) = Uo, 1 < x < N x /4, 3N x /4 < x < N x , 

u y (x, 0) = —Uo, N x /4 < x < 3N x /4. 

The initial condition u x (x, 0) is set to zero every where. We consider two separate cases with and without 
a constant streaming velocity V . 

7.2. Steady case ( V = 0). For the case of zero streaming velocity, the initial condition for u x is zero 
in the system. The solution of the Navier-Stokes equation for this simple problem is: 

(7.6) u y (x, t) = a n exp {-v n k^t) cos (k n x) , 
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X X 

Fig. 4. Decay of discontinuous shear wave velocity profile u y (x, t). The lines and symbols (X ) are theoretical [Eq. (7.6)] 
and numerical results, respectively. Only the positive half of the velocity profiles are shown in the Figures, (a) The LBE model 
with no interpolation, (b) with the central interpolation and r = 0.5. 


where a n are the Fourier coefficients of the initial velocity profile u y (x, 0), v n = u(k n ), and k n — 2ir(2 n — 
1 )/N x . The magnitude of the u y (x, 0), Uo = 0.0001 in the simulations. 

Figures 4(a) and 4(b) show the decay of the rectangular shear wave simulated by the normal LBE scheme 
and the LBE scheme with second-order central interpolation (with r = 0.5, where r is the ratio between 
advection length 6 X and grid size A^), respectively. (The detailed analysis of LBE schemes with various 
interpolations is provided in Appendix B.) The lines are theoretical results of Eq. (7.6) with v(k n ) obtained 
numerically. The times at which the profile of u y (x, t ) (normalized by Uo) shown in Figs. 4 are t = 100, 200, 
. . . , 500. The numerical and theoretical results agree closely with each other. The close agreement shows the 
accuracy of the theory. In Fig. 4(b), the overshoots at early times due to the discontinuous initial condition 
are well captured by the analysis. This overshoot is entirely due to the strong fc-dependence of v(k) caused 
by the interpolation. This phenomena is not necessarily connected to the Burnett effect, as claimed by a 
previous work [38]. This artifact is also commonly observed in other CFD methods involving interpolations. 

Figure 5 shows the decay of u y (x, t ) at one location of discontinuity, x = 3N x /4 = 63. We tested the 
normal LBE scheme without interpolation and the LBE scheme with second order central interpolation with 
r = 0.5, and compared the numerical results with theoretical ones. Again, the numerical and theoretical 
results very well agree with each other for both cases (with and without interpolation). Note that the time 
is rescaled as r~ 2 t in the Figure. It should be pointed out that the LBE solutions of the flow differ from 
the analytic solution of the Navier-Stokes equation in both short time and long time behavior. Interpolation 
causes overshoot in the velocity at the initial stage. Even without interpolation, the LBE solution does not 
decay (exponentially) right away. This is due to the variation of the viscosity with k and this could be 
interpreted as the influence of the kinetic modes. (If we had a vanishingly small Knudsen number then the 
fc-dependence would be negligible, however all relaxation rates must be smaller than 2 so that higher modes 
can play a role). This transient behavior is due to the higher order effect (of velocity gradient), as discussed 
previously. 

7.3. Streaming case ( V = constant). We also consider the case with a constant streaming in the 
initial velocity, i.e., u x (x, 0) = V x = 0.08. This allows us to check the effects of the non-Galilean invariance 
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Fig. 5. Decay of discontinuous shear wave velocity u y (x, t) at a location close to the discontinuity x = ZN X /A. The solid 
lines and symbols (x) are theoretical and numerical results, respectively. The LBE scheme with no interpolation does not have 
a overshooting, whereas the LBE scheme with the central interpolation and r = 0.5 has. The time is rescaled as r~“t. 


in the system. With a constant streaming velocity, the solution of the Navier-Stokes equation is: 

(7.7) u y (x, t ) = ^ a n exp (-v n k 2 n t) cos[k n (x - g n V x t )] , 

n 

where g n = g(k n ) is the Galilean-coefficient. 

Similar to Figs. 4, Figs. 6 show the evolution of u y {x, t) for the same times as in Figs. 4. The solid lines 
and the symbols ( X ) represent theoretical and numerical results, respectively. Shocks move from left to right 
with a constant velocity V x = 0.08. Figs. 6(a), 6(b), and 6(c) show the results for the normal LBE scheme 
without interpolation, the scheme with second order central interpolation, and the scheme with second order 
upwind interpolation, respectively. In Figs. 6(b) and 6(c), the dotted-lines are the results obtained by setting 
g n = 1 in Eq. (7.7). Clearly, the effect of g(k) is significant. For the LBE scheme with central interpolation, 
the results in Fig. 6(b) with g(k ) = 1 under-predict the overshooting at the leading edge of the shock and 
over-predict the overshooting at the trailing edge; whereas the results in Fig. 6(c) for the LBE scheme with 
upwind interpolation over-predict the overshooting at the leading edge of the shock and under-predict the 
overshooting at the trailing edge. 

8. Conclusion and Discussion. In this paper, a generalized 9-velocity LBE model based on the 
generalized LBE model of d’Humieres [8] is presented. The model has the maximum number of adjustable 
parameters allowed by the discrete velocity set. The value of the adjustable parameters are obtained by 
optimizing the hydrodynamic properties of the model through the linear analysis of the LBE evolution 
operator. The linear analysis also provides the generalized hydrodynamics of the LBE model, from which 
dispersion, dissipation, isotropy, and stability of the model can be easily analyzed. In summary, a systematic 
and general procedure to analyze the LBE models is described in detail in this paper. Although the model 
studied in this paper is relatively simple, the proposed procedure can be readily applied to analyze more 
complicated LBE models. 

The theoretical analysis of the model is verified through numerical simulation of various flows. The 
theoretical results closely predict the numerical results. The stability of the model is also analyzed and 
compared with the BGK LBE model. It is found that the mechanism of separate relaxations for the kinetic 
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Fig. 6. Decay of discontinuous shear wave velocity profile u y (x , #) with a constant streaming velocity V x = 0.08. The 
thick lines and symbols (x) are theoretical [Eq. (7.7)] and numerical results, respectively. The thin lines in (b) and (c) are 
obtained by setting g n = 1 in Eq. (7.7). (a) The LBE model with no interpolation, (b) with the central interpolation and 
r = 0.5, (c) with the upwind interpolation and r = 0.5. 

modes leads to a model which is much more stable than the BGK LBE model. 

The proposed model is a Galerkin type of scheme. In comparison with the BGK LBE model, the 
proposed model requires the transformations between the discrete velocity space V and the moment space 
M back and forth in each step in the evolution equation. However, the extra computational cost due to 
this transformation is only about 10 - 20% of the total computing time. Thus, the computational efficiency 
is comparable to the BGK LBE model. Our analysis also shows that the LBE models with interpolation 
schemes have enormous numerical hyper-viscosities and anisotropies due to the interpolations. 

We also find optimal features of the proposed 9-velocity model: it is difficult to improve the model 
by simply adding more velocities. For instance, we found that adding eight more velocities (±1, ±2) and 
(±2, ±1) would not improve the isotropy of the model. However, our analysis does not provide any a priori 
knowledge of an optimal set of discrete velocities. That problem can only be solved by optimization of the 
moment problem in velocity space [20]. It is also worth noting that the values of the adjustable parameters 
in our model coincide with the corresponding parameters in the BGK LBE model except two (c*3 and 74). 
The main distinction between our model and the BGK LBE model is that in our model has the freedom 
to allow the kinetic modes to relax differently, whereas in the BGK LBE model, all kinetic modes relax at 
the same rate. This mechanism severely affects the stability of the BGK LBE schemes especially when the 
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system is strongly over-relaxed. 

It should be mentioned that the procedure we propose here can be applied to analyze the linear stability 
of spatially nonuniform flows, such as Couette flow, Poiseuille flow, or lid-driven cavity flow. For spatially 
nonuniform flows, the lattice Boltzmann equation is linearized over a finite domain including boundary 
conditions. This leads to an eigenvalue problem with many more degrees of freedom as was needed in the 
analysis of this paper. Standard Arnoldi techniques [39] allow us to determine parts of the spectrum of the 
linearized collision operator, in particular to study the flow stability. This analysis enables us to understand 
the observation that some flows are much more stable than what is predicted by the linear analysis of 
spatially uniform flows. For instance, in plane Couette flow with only 2 nodes along the flow direction, the 
only possible values of k along the same direction are 0 and n, which are far away from the value of k at 
which the bulk instability occurs. Namely, the reciprocal lattice fc is not large enough to accommodate the 
possible unstable modes. Furthermore, in the direction perpendicular to the flow, although the reciprocal 
lattice k can accommodate unstable shear modes, the velocity gradient, however, alters the stability of the 
system. (It improves the stability in this particular case.) 

One philosophic point must be stressed. We deliberately did not derive the macroscopic equations 
corresponding to the LBE model in this work; instead, we only analyzed the generalized hydrodynamic 
behavior of the modes of the linearized LBE evolution operator. We argue that if the hydrodynamic modes 
behave exactly the same way as those of the linearized Navier-Stokes equations, up to a certain order of fc, 
provided that the Galilean invariance is also assured up to a certain order of fc, then we can claim that the 
LBE model is indeed adequate to simulate the Navier-Stokes equations (up to a certain order of fc). There 
is no distinction between the LBE model and the Navier-Stokes equations up to a certain order of fc. Thus, 
there is no need to use the Chapman-Enskog analysis to obtain the macroscopic equations from the LBE 
models. On the other hand, we have also shown that, in the limit of fc = 0 , these two approaches obtain 
the same results in terms of the transport coefficients and the Galilean-coefficient. Nevertheless, it is very 
difficult to apply the Chapman-Enskog analysis to obtain the generalized hydrodynamics of the LBE models, 
which is important to LBE numerical simulations of hydrodynamic systems. The stability result obtained 
by the linear analysis presented in this paper is very difficult for the standard Chapman-Enskog analysis to 
obtain. Therefore, the proposed procedure to analyze the LBE model indeed contains more information and 
is more general than the low order Chapman-Enskog analysis. Albeit its generality and powerfulness, the 
linear analysis has its limitations. Because it is a local analysis, it does not deal with gradients. 

Our future work is to extend the analysis to fully thermal and compressible LBE models in three- 
dimensional space. 
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Appendix A. Coupling Between Density and Other Modes. 

To consider the coupling between the density fluctuation 5p = p — (p) and other modes, e, s, p xx , and 
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p xy , the equilibrium values of these modes are modified as the following: 


(A. la) 

e (eq) = a 2 P + 72 (i 2 + jy)( 2 - p), 

(A.lb) 

£ (eq) = c* 3 p + 74 (j 2 + jy)( 2 - p), 

(A.lc) 

Pxx ] = 7i til + jy)(2 ~ p), 

(A.ld) 

Px : = 73 tixjy){ 2 - p). 


where (2 — p) is used to linearly approximate 1/p when the averaged density po = (p) = 1. With the above 
modifications, four elements in the first column of the linearized collision operator C accordingly become: 


(A. 2a) 
(A. 2b) 
(A. 2c) 
(A. 2d) 


C 

c 

c 

c 


12 — S 2 


13 = S 3 


\c*-\'t2<y? + v;) 
4 a 3 - + v f) 


lS = ~\ssll{V 2 -V 2 ), 

19 = — 2 S 9T3^4^3 i- 


Based on the linearized collision operator with the above changes, the shear and the bulk viscosity at 
the limit of k ->• 0 are: 


(A.3) v 0 = ^(1 - V 2 cos 2 <j>) ^ ^ , 

Vo = T ^2 ( 2 " 3c «X 2 - * 2 ) ~ (! - 3c 2 )(3s 2 s 8 - 2 s 2 - 4s 8 ) 

\2tS‘2 l^CgS2^8 

V 2 

+ - [s 2 - S S + 2(s 2 .S 8 - S 2 ~ Sg) COS 2 (/>] 

4s 2 s 8 

, . V s cos 6, , 9 ,, 

(A. 4) +- — [s 2 - sg + .5 2 (s 8 - 2) cos 2 <f>\ . 

4 c s $ 2 Sg 

The sound modes propagate with velocity V ± c s (at first order in k). The Galilean-coefficient up to 0(k 2 ) 
is: 


(A.5) 




[(s 8 - 2)(s 5 - S 8 )(S5«8 - 3s 5 - 3s 8 + 6) + (cos 4 9 - cos 2 9)] 
[(2 - s 8 )(sg - s 2 ) sin 2 + 2c 2 s 2 (s| - 6s g + 6) cos 2 <j>]. 


Appendix B. Interpolated LBE Scheme. 

Recently, it has been proposed to use interpolation schemes to interpolate {/„} from a fine mesh to a 
coarse mesh in order to improve the spatial resolution calculations for a limited cost in total number of nodes 
[21, 22]. Obviously, the interpolation schemes create additional numerical viscosities. The Chapman-Enskog 
analysis shows that any second or higher order interpolation scheme does not affect the viscosities in the 
limit k — ¥ 0 on the fine mesh. A problem with much greater importance in practice is to calculate the 
viscosity at finite k. To our knowledge, no such analysis is now available in the literature. 

In the interpolated LBE schemes, the advection step is altered by the interpolation scheme chosen, 
whilst the collision step remains unchanged. The advection on a fine mesh combined with interpolation on 
a coarse mesh is the reconstruction step on the coarse mesh. Therefore, to obtain the modified linearized 
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evolution operator L only the advection operation A must be changed. In what follows, we shall consider 
a coarse mesh with lattice constant S x , and time step St ■ The lattice constant of a underlying fine mesh is 
rS x , with r < 1. Effectively, the hopping velocities of particles are reduced by a factor of r on coarse mesh. 
Therefore, dimensional analysis suggests that the sound speed is reduced by a factor of r, and the viscosities 
are reduced by a factor of r 2 in the limit k = 0. However, the dimensional analysis does not provide any 
information about the quantitative effects of interpolation when k is finite. We shall analyze the effects of 
some commonly used second-order interpolation schemes in the LBE methods. For simplicity, we shall only 
deal with a uniform mesh with square grids. 

B.l. Central interpolation. The reconstruction step with second-order central interpolation is given 
by the following formula: 

(B.l) f a (r j ) = r -^p±f*( rj - Sr a ) + (1 - r 2 )f*(r j) + 1 ^±Vf*(r j + Sr a ) , 

where /* is the post-collision value of f a , i.e., 

(B-2) f* = f Q + njf a ) , 

and 

(B.3) Sr a = - e a . 

r 

The advection operator in this case becomes: 


(B.4) 


A = diag(l, A, C, B, D, AC, CB, BD, DA) , 


where 

(B.5a) 

(B.5b) 

(B.5c) 

(B.5d) 


A = 
B = 
C = 
D = 


r(V + 1)P I fl r 2 l I r(r ~ 1} 

2 r ,+ 2p , 

r ( r + , (1 _ r 2x , r(r - 1 ) p 

2p >+ 2 ’ 

r(r + l)q . 2 , r(r - 1) 

2 r ,+ 2g . 

r ( r + !) . n _ 

2 q >+ 2 ’ 


where p = e lkx and q = e lky . With the new phase factors, we find new results at order 1 and 2 in k. The 
speed of sound and the “Galilean coefficient” are multiplied by r and the viscosity coefficients are multiplied 
by r 2 . 

At higher order in k, dispersion effects due to lattice arise, leading to differences between solutions of 
the standard Navier-Stokes equations and the flows computed used the LBE technique. 

As in Eq. (5.24), we find that the advection coefficient for shear waves can be made isotropic to second 
order in k by choosing 


(B.6) 


s 5 = 3r 2 


(2 - s 8 ) 

(3 r 2 - s s ) ’ 


which improves Eq. (5.24), since we can choose Sg close to 2 while maintaining s 5 reasonably far away from 
2 (between 1 and 3/2) by taking r 2 close to 2/3. 
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B.2. Upwind interpolation. The upwind direction in the LBE method is relative to the particle 
velocity e a (the characteristics) rather than the flow velocity u. Therefore, the interpolation stencil is static 
in time. Second-order upwind interpolation leads to 

(B.7) f a (vj) = r ±^±f*( rj -2 6r a ) + r(2 - r)/>, - 5r a ) + (1 ~ r )( 2 ~ r ) f * {r .) , 

where 5r a is defined in Eq. (B.3). Accordingly, the phase factors in the advection operator given by Eq. (B.4) 
become: 


(B.8a) 

(B.8b) 

(B.8c) 

(B.8d) 


^ (l-r)(2-r) , r(2 — r) , r(r - 1) 

2 + p + 2p2 ’ 

B= (l-rK2-r) +r(2 _ r)i)+ r(r-lV ^ 
^ (1— r)(2 — r) , r(2-r) , r(r - 1) 

z)= (l-rK2-r) +r(2 _ r)g+ r(r-ltf 


where p = e lkx and q = e tky . 

Again, the third order term ( g \ ) in k for the shear mode is anisotropic unless the following relation is 
satisfied: 


(B-9) 


.2 '^8 ) 

(3r 2 -3rs 8 + 2s 8 )' 


For sg and S 5 in the usual range (s§ near 2 and s 5 between 1 and 3/2), the preceding equation leads to a 
complex value of r. It should be pointed out that due to the commutativity of propagation along x-axis and 
y-axis, one could apply different interpolation formulae along each axis, according to physics of the flow. For 
instance, large stretch of grid can be applied in the direction along which flow fields do not change much 
in space, whereas in the other orthogonal direction, a. normal grid (without interpolation) or even a refined 
grid [9] can be used, so that the aspect ratio of the meshes is large enough to be appropriate to the flow. 

Figures 7 show the fc-dependence of the normalized shear viscosity v{k ) /vq for the LBE model with and 
without interpolation schemes. Three orientations of k are chosen: 0 = 0 (solid line), 7 t/ 8 (dotted line), and 
7r/4 (dashed line). Figs. 7(a), 7(b), and 7(c) show the v(k)/vo for the LBE model with no interpolation, with 
second order central interpolation scheme and r = 0.5, and with second order upwind interpolation scheme 
and r = 0.5, respectively. It should be stressed that interpolation schemes do create an enormous amount 
of numerical viscosity at k = 7t/ 2: Both the central and the upwind interpolation schemes increase the 
shear viscosity at k = 7t/2 by almost two order of magnitude, whereas without interpolation, corresponding 
increase for the LBE scheme is at most only a factor of about 2.5 (in the direction 0 = 7t/ 8). In all cases, 
the viscosity displays significant anisotropy at k = n/ 2 . 

Similar to Figs. 7, Fig. 8 shows the fc-dependence of the Galilean-coefficient g(k). The three curves in 
the middle of the figure corresponding to the LBE model without interpolation. The lower three curves, 
g(k) < 1, correspond to the LBE scheme with the central interpolation, and the upper three curves, g(k) > 1, 
correspond to the LBE scheme with the upwind interpolation. Again, interpolations create significant effects 
on Galilean invariance. 


One common feature observed in Figs. 7 and 8 is that the transport coefficients of a model along the 
direction of 0 = ir /8 is far apart from those along the directions 9 = 0 and 9 = 7t/4. This is related to the 
fact that for the square lattice, the wave vector k along the direction 9 = 7r/8 is not a reciprocal vector of 
the underlying lattice. 
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0 rr/4 tt/2 

k 

Fig. 7. The k-dependence of viscosities for various models. The values of the adjustable parameters and the relaxation 
parameters are the same as in Figs. 1. The solid lines, dotted lines, and dashed lines correspond to 0 = 0, tt/8, and tt/4, 
respectively, (a) The normal LBE model with no interpolation, (b) with central interpolation, and (c) with upwind interpolation. 



Fig. 8. The k-dependence of the Galilean-coefficient g for various models. The solid lines, dotted lines, and dashed lines 
correspond to 6 = 0, n/8, and n/4, respectively. The middle three curves are g(k) for the LBE model without interpolation, 
the lower three for the LBE model with the central interpolation and r = 0.5, and the upper three for the LBE model with the 
upwind interpolation and r = 0.5. 
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